function y = noiseReduction(x, FS)

    %Fecha:          01/09/08
    %MATLAB R2008a UNIX-Linux 64

    %Syntax:        cleanSpeech = noisereduction(x, FS);
    %FS es frecuencia de muestreo

    %Esta función acepta la señal de voz ruidosa y su frecuencia de muestreo como
    %entradas y devuelve la versión 'denoised' de la señal.

    %Basada en Sustraccion Espectral de Boll (1979).
    %Referencia:
    %S. F. Boll, “Suppression of Acoustic Noise in Speech using
    %Spectral Subtraction,” IEEE Trans. on ASSP, vol. 27(2), pp.
    %113–120, 1979.


speech = x;                      %Recorta ruido y artefactos. Llama a la funcion Voice Activity Detection.



%Definicion de constantes

le          = length(speech);                                               %Total de muestras de la señal de voz
s_time      = 0.1;                                                          %Un periodo de silencio inicial en seg
w_time      = 0.025;                                                        %Tamaño de los cuadros en segundos.
olap_time   = 0.005;                                                        %Solapamiento en segundos
w           = fix(w_time*FS);                                               %No de Muestras/ventana (Samples/win) 
olap        = fix(olap_time*FS);                %No de muestras/solapamiento
s           = fix((s_time-w_time)/olap_time)+1; %Número de segmentos de silencio 
tsg         = fix((le-w)/olap)+1;               %Número de segmentos en la señal. tsg 
win         = hann(w).^ 0.5;
%win        = hamming(w);                       %Ventana Hamming
%win        = hann(w);                          %Ventana Hanning
%win        = blackman(w);                      %Ventana Blackman
%win        = bartlett(w);                      %Ventana Bartlett
%win        = gausswin(w);                      %Ventana gauss, alfa= 2,5
xofT        = zeros(tsg,w);                     %Matriz para mantener los segmentos


%Segmentación de una sección de voz (speech)
for     i = 1:tsg
            nextIndex = (i-1)*olap+1;
            xofT(i,:) = speech(nextIndex:nextIndex+w-1).*win;
end


%Matriz de segmentos de silencios
s_seg = zeros(s,w);
for     k = 1:s
            nextIndex = (k-1)*olap+1;
            s_seg(k,:) = x(nextIndex:nextIndex+w-1).*hamming(w);
end


%Se obtienen los espectros del habla y ruido, también la fase de la señal
xofK    = fft(xofT,w,2);
halfn   = 1+floor(w/2);
xofK    = xofK(:,1:halfn);
Phase   = angle(xofK);
nofK    = fft(s_seg,w,2);
nofK    = nofK(:,1:halfn);
nmean   = mean(nofK);                           %Calcula la media del espectro de ruido


%Promedia los valores de 3 cuadros sucesivos del espectro de voz
xofK1 = xofK;
for     l = 2:tsg-1
        xofK1(l,:) = (xofK(l-1,:)+xofK(l,:)+xofK(l+1,:))/3;
end


[nr,nc] = size(xofK1);
nmean   = repmat(nmean,nr,1);
xofK1   = xofK1-nmean;               %Sustraccion espectral (SS)
xofK2   = max(xofK1,zeros(nr,nc));   %Rectificación de media onda. Remueve valores -ve >> valores negativos

xoKwP   = xofK2.*exp(j*Phase);       %Multiplica estimación de voz con fase ¡?
xoKwPc  = conj(xoKwP);


%Dobla los espectros conjugados alrededor de los espectros principales para formar una FFT completa
if      isinteger(nc/2)
            xoKwP = [xoKwP,fliplr(xoKwPc(:,2:end-1))];
else
            xoKwP = [xoKwP,fliplr(xoKwPc(:,2:end))];
end


denoisedSignal = zeros((nr-1)*olap+w,1);%En esta variable se mantendrá el resultado final


%Obtiene la IFFT real y hace adición de los solapamientos, luego de ventanear con 'win'
for 	m = 1:nr
	nextIndex = (m-1)*olap+1;
	invfft    = real(ifft(xoKwP(m,:),w))'.*win;
	denoisedSignal(nextIndex:nextIndex+w-1) = denoisedSignal(nextIndex:nextIndex+w-1)+invfft;
end

y=myVAD(denoisedSignal);

% Representacion Grafica. Front End de Reduccion de Ruido y Recorte de
% Silencios

% totalduration   =   length(denoisedSignal)/FS;
% timeduration    =   0:(totalduration/length(denoisedSignal)):totalduration-(totalduration/length(denoisedSignal));
% totalduration2  =   le/FS;
% timeduration2   =   0:(totalduration2/le):totalduration2-(totalduration2/le);
% 
% 
% 
% subplot(2,2,3);
% 	plot(timeduration,denoisedSignal); title('Señal con reduccion de ruido', 'FontSize', 12);
% subplot(2,2,4); 
% 	spectrogram(denoisedSignal); title('Espectrograma de la Señal con reducción de ruido', 'FontSize',12);
% subplot(2,2,1)
% 	plot(timeduration2,speech); title('Señal de Entrada', 'FontSize', 12);
% subplot(2,2,2)
% 	spectrogram(speech); title('Espectrograma de la Señal de entrada', 'FontSize', 12);
    

    
    
    
    

